## medina (2018) or fieler (2011) style utility function
##	non-homothetic, CES

## continuum of multiproduct firms


require(data.table)
require(ggplot2)

theme_dan <- function(base_size=12){
  	theme_classic(base_size = base_size)+
	theme(axis.line.x = element_line(color='black'), axis.line.y=element_line(color='black')) + 
	theme(legend.position='bottom', legend.title=element_blank())  
}


## =========================================================================
## PROFIT COMPUTATION
## =========================================================================

lambda <- function(x, a, sig, p, Y, N, J){
	return(sum((x^(-sig))*a*N*J*(p^(1 - sig))) - Y)
}

#uniroot(lambda, interval=c(0.1, 5), a=a, sig=sig, p=c(3,2), Y=Y, N=N)$root
#plot(sapply(seq(1, 2, length.out=20), function(x) lambda(x, a, sig, p=c(3, 2), Y=10, N=c(10,10))))

profit <- function(a, sig, mc, N, J, g, Y, lambda.f){
	## input:
	##	a: quality shifters
	##	sig: CES parameters
	##	mc: marginal costs
	##	N: number of each type of SKU
	##	g: ER shifter
	##	Y: income
	## output:
	##	profit for each type of SKU

	p <- g*mc*(sig/(sig - 1))		## individual price optimizations
	la <- uniroot(lambda.f, interval=c(1e-06, 1000), a=a, sig=sig, p=p, Y=Y, N=N, J=J)$root
	q <- (la^(-sig))*(p^(-sig))*a
	pi <- q*(p - mc*g)

	return(pi)
}

gseq <- seq(1,3,length.out=10)
piout <- sapply(1:length(gseq), function(x) profit(a=a, sig=sig, mc=mc, N=c(2,2), J=J, g=gseq[x], Y=Y, lambda.f=lambda))


## =========================================================================
## ENTRY MODEL
## =========================================================================

## firms have quadratic costs in the number of each kind of product (h and l)
## maximize N_m*pi_m - 0.5*f*N_m^2
##	=> N_m = pi_m/f

Nsim <- function(a, sig, mc, f, g, Y, J, lambda.f){
	## compute profit (for h and l type products) for each firm, solve
	## 	for optimal number of products

	N0 <- c(2,2)
	N1 <- N0 + 1
	
	iter <- 0
	while(max(abs(N0 - N1)) > 0.01 & iter < 100){
		iter <- iter + 1
		N0 <- N1
		pivar <- profit(a, sig, mc, N=N1, J, g, Y, lambda.f=lambda.f)
		N1[1] <- ifelse(pivar[1]/f > N0[1], min(pivar[1]/f, N0[1] + 0.1), max(pivar[1]/f,N0[1] - 0.1))
		N1[2] <- ifelse(pivar[2]/f > N0[2], min(pivar[2]/f, N0[2] + 0.1), max(pivar[2]/f,N0[2] - 0.1))
	}
	return(N1)
}

## =========================================================================
## SIMULATION
## =========================================================================

## SET OUTPUT DIRECTORY HERE 
## setwd(...)

## parameters
Y <- 10
a <- c(7.5, 2.5)
sig <- c(3, 2.5)
mc <- c(2.75, 2.5)
f <- 0.1
J <- 10			## number of firms

gseq <- seq(1,2,length.out=23)
Nout <- sapply(1:length(gseq), function(x) Nsim(a=a, sig=sig, mc=mc, f=0.1, g=gseq[x], Y=Y, J=J, lambda.f=lambda))
piout <- sapply(1:length(gseq), function(x) profit(a=a, sig=sig, N=Nout[,x], mc=mc, g=gseq[x], Y=Y, J=J, lambda.f=lambda))


simdat <- data.table(er = c(gseq, gseq), quality = c(rep('High Quality', length(gseq)), rep('Low Quality', length(gseq))), 
						pi = c(piout[1,], piout[2,]), N = c(Nout[1,], Nout[2,]))

p1 <- ggplot(simdat, aes(x = er, y = pi, color = quality, group = quality, linetype=quality)) + geom_smooth(se=F) + 
	xlab('Normalized exchange rate') + ylab('Profit') + theme_dan()

p2 <- ggplot(simdat, aes(x = er, y = N, color = quality, group = quality, linetype=quality)) + geom_smooth(se=F)+ 
	xlab('Normalized exchange rate') + ylab('Number (mass) of products') + theme_dan()


pdf('fieler_utility_profit_simulation.pdf', height=4, width=4)
p1
dev.off()

pdf('fieler_utility_product_choice_simulation.pdf', height=4, width=4)
p2
dev.off()





































##